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its rate of dissipation (mVsVs (ft /s“/s)), 
respectively 

K an acceleration parameter 

k gas thermal conductivity, W/m-K (Btu/hr-ft- 

R) 

/ mixing length constant usually known as 

VonKarman length scale, m (ft) 

m mass flow rate, kg/s (lbm/s) 

MW molecular weight, kg/kgmole 

M Mach number 

Md mass fraction of particles with diameter 

greater than D 

N particles spread diameter 

P 0t P total and static pressure, respectively, Pa 

(psia) 

Pr Prandtl number 

Pr t turbulent Prandtl number 

non-dimensional pressure gradient 

q wall heat flux, W/rrf (Btu/hr-ft 2 ) 

Re t i Reynolds number, p d/ ju 

9i recovery factor 

r p propellant bum rate, m/s (ft/s) 

r radius, m (in.) 

S source term (Eq. 1 ) 

T 0t T total and static temperature, respectively, K 

(R) 

non-dimensional temperature in wall 
coordinates 

u, v axial and radial velocity components, m/s 
(ft/s) 

tt M velocity magnitude of a and v at the motor 


centerline (m/s, ft/s) 

if non-dimensional velocity in wall 
coordinates 

u ' velocity fluctuation, m/s (ft/s) 

u T shear velocity, m/s (ft/s) 

x , y axial and radial axes, m (in.) 

y + non-dimensional distance from the wall in 
wall units 

y ratio of specific heats, Q/C 

p , v dynamic (N-s/m 2 (Ibm/ft-s)) and kinematic 

(m 2 /s (ft 2 /s)) viscosity, respectively 
p gas density, kg/m 3 (lbn/ft 3 ) 

0 Inlet flow direction (0 = 0 degree is normal) 
Xw wall shear stress (N/m" (lbf/ft 2 )) 

(y) function of radial direction y 

Subscripts 

a rate of accretion flux 

roo recovery temperature 

c concentration 

DP discrete phase 

e at nozzle exit 

g gas 

h hydraulic 

o chamber conditions 

p particles, propellant 

ref reference 

s static, surface 

vac vacuum 

w nozzle wall 

oo conditions at motor centerline 

1 , 2 based on Method 1 or Method 2, also based 
on internal and external boundary layer. 


o 
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Abstract 

An enhanced performance solid rocket booster concept for the space shuttle system has been proposed. The 
concept booster will have strong commonality with the existing, proven, reliable four-segment Space Shuttle Reusable 
Solid Rocket Motors (RSRM) with individual component design (nozzle, insulator, etc.) optimized for a five-segment 
configuration. Increased performance is desirable to further enhance safety/reliability and/or increase payload 
capability. Performance increase will be achieved by adding a fifth propellant segment to the current four-segment 
booster and opening the throat to accommodate the increased mass flow while maintaining current pressure levels. 
One development concept under consideration is the static test of a "standard” RSRM with a fifth propellant segment 
inserted and appropriate minimum motor modifications. Feasibility studies are being conducted to assess the 
potential for any significant departure in component performance/loading from the well-characterized RSRM. An 
area of concern is the aft motor (submerged nozzle inlet, aft dome, etc.) where the altered internal flow resulting from 
the performance enhancing features (25% increase in mass flow rate, higher Mach numbers, modified subsonic nozzle 
contour) may result in increased component erosion and char. To assess this issue and to define the minimum design 
changes required to successfully static test a fifth segment RSRM engineering test motor, internal flow studies have 
been initiated. Internal aero-thermal environments were quantified in terms of conventional convective heating and 
discrete phase alumina particle impact/concentration and accretion calculations via Computational Fluid Dynamics 
(CFD) simulation. Two sets of comparative CFD simulations of the RSRM and the five-segment (FSM) concept motor 
were conducted with CFD commercial code FLUENT. The first simulation involved a two-dimensional axi-symmetric 
model of the full motor, initial grain RSRM. The second set of analyses included three-dimensional models of the 
RSRM and FSM aft motors with four-degree vectored nozzles. 

d diameter, m (in.) 


Nomenclature 

A area, m" (ft") 

A* Van Driest constant 

C f friction coefficient 

Cp,Cv specific heats at constant pressure and 

volume, respectively, J/kg-K (Btu/lb m -R) 


F thrust, N (lbf) 

g acceleration due to gravity, m/s 2 (ft/s 2 ) 
h convective heat transfer coefficient, W/m 2 - 

K (Btu/hr-fr-R) 

/ turbulence intensity, a'/ii (%) 

l sp specific impulse, s (lbf-s/lbm) 

K, e turbulence kinetic energy (m 2 /s 2 (ft7s 2 )) and 
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respectively 

Superscripts 

n bum rate pressure exponent 

* throat conditions 

flux 

— mean and time averaged 

Introduction: 

Internal flow simulations of the proposed five 
segments motor were conducted in order to assess 
the aft motor aero-thermal environment differences 
between the five-Segment Motor (FSM) and the 
four-Segment Redesigned Solid Rocket Motor 
(RSRM). Design differences between the two 
systems indicated that any significant departure in 
component erosion/char from the RSRM would 
likely be manifest in the submerged and nozzle 
regions. Objectives of these analyses are: 

• Compare the predicted aero-thermal 
environments in the aft motor region and 
nozzle. 

• Identify any potential “show stoppers” 
associated with the five-segment concept. 

• Flag those areas or phenomena that should be 
further scrutinized in the development phase 
of the program. 

The five segments booster was proposed as an 
enhanced performance solid rocket booster concept 
for the space shuttle system. The concept booster 
will have strong commonality with the existing, 
proven, reliable four-segment RSRM's with 


individual component design (nozzle, insulator, etc.) 
optimized for a five-segment configuration. 
Increased performance is desirable to further 
enhance safety/reliability and/or increase payload 
capability. Performance increase will be achieved 
by adding a fifth propellant segment to the current 
four-segment booster and opening the throat to 
accommodate the increased mass flow while 
maintaining current pressure levels. 

Feasibility studies are being conducted to assess 
the potential for any significant departure in 
component performance/loading from the well- 
characterized RSRM. An area of concern is the aft 
motor (submerged nozzle inlet, aft dome, etc.), 
where the altered internal flow resulting from the 
performance enhancing features (25% increase in 
mass flow rate, higher Mach numbers, modified 
subsonic nozzle contour) may result in increased 
component erosion and char. Internal aero-thermal 
environments were quantified in terms of 
conventional convective heating and discrete phase 
alumina particle concentration and accretion 
calculations via CFD simulation. 

Discussion: 

The following literature survey is a discussion of 
what has and has not been dealt with in the world of 
solid rocket motors. Selected studies are discussed 
in the following paragraph and are referenced 
chronologically in [1-29]. Slag deposition process 
through flow modeling and the subsequent 
accumulation and pooling of slag material within the 
Space Shuttle solid rocket motor (SRM) was first 



documented by Boraas [1]. Flow field modeling 
consisted of a potential flow gaseous computation 
with a Lagrangian particle trajectory calculation. 
The separated flow region in the aft end of the 
submerged nozzle motor was approximated by a 
fictitious interface drawn between the submerged 
nozzle entrance and the aft dome recirculation 
region. Particle trajectories were computed by 
varying particle starting locations on the propellant 
surface. Murdock [2] discussed the ejection of a 
body that blocks the nozzle throat and may cause 
perturbation in the measured thrust-time histories of 
solid rocket motors. Traineau at al. [3] conducted 
cold-flow simulation tests of a nozzleless solid 
rocket motor using a two-dimensional porous-walled 
duct with an impermeable-walled diverging section. 
Beddini [4] conducted theoretical analysis of the 
flow in porous-walled tubes and channels with 
appreciable injection through the duct wall. The 
flow at large injection Reynolds numbers can 
undergo at least three flow regimes. Basset [5] 
compiled results of modeling of circumferential 
flow in the RSRM induced by potential non- 
axisymmetric flow sources. It describes analytical 
and experimental results. The numerical techniques 
utilized eight 3D fluid-flow codes and two 3D heat 
flow codes, while the testing media included water, 
cold air, and hot combustion gas. Salita [6] 
suggested a bimodal log-normal droplet size 
distribution at a pressure of about 6.90 MPa (1000 
psia) consisting of small droplets with 1.5 pm mean 
and large droplets with a 100 pm mean. Whitesides 
et al. [7] conducted a series of subscale cold flow 


tests to quantify the gas flow characteristics at the 
aft end of the RSRM. In specific, measurements of 
static pressure and gas velocities in the vicinity of 
the nozzle/case joint when nozzle is gimbaled at 
angles of 0, 3.5, and 7 degrees. Golafshani and Loh 
[8] conducted a time-dependent, axisymmetric 
numerical solution of the Navier-Stokes to analyze 
the viscous coupled gas-particle non-reacting flow 
in solid rocket motors. The solution assumed 
laminar internal flow. Waesche et al. [9] conducted 
flow-visualization tests in a 1 /8-scale model of the 
RSRM in a water tunnel to simulate (1) 
circumferential flow induced by asymmetric 
inhibitor stub in the port of the RSRM, (2) vortex 
formation in the aft-dome cavity with and without 
nozzle vectoring. Burning propellant was simulated 
through the introduction of a uniform distribution of 
water along simulated bum-back patterns. Vortices 
shed from protruding inhibitor were found to 
diminish by wall injection. Circumferential flow 
resulting from a flawed inhibitor was limited to 
region near the missing portion of the inhibitor. 
Strong circumferential flow in the aft dome was 
observed toward the end of bum. Boraas [10] used 
a water table to simulate flow in the aft dome of the 
RSRM. A thin layer of water was pumped across a 
table between sidewalls that simulated the motor 
centerline and its aft-dome-nozzle boundary in non- 
vectoring condition and late in bum. Lightweight 
dye was injected uniformly across the flow field at 
an upstream location. Entrainment of the dye into 
the cavity revealed a counter-rotating three vortices 
in the aft dome. Majumdar et al. [11] conducted 
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three-dimensional CFD analyses to calculate 
circumferential pressure and velocity gradients in 
the vicinity of an asymmetric inhibitor stub in the 
port of the RSRM. The numerical predictions were 
compared with the measurements from a 7.5% scale, 
cold-flow model of the RSRM. A maximum of 0.34 
atm (5 psid) has been shown. Waesche et al. [12] 
conducted a series of cold-gas tests in a 1 /8-scale 
(18-in. dia.) model of the forward segment of 
candidate grain designs for the Advanced Solid 
Rocket Motor (ASRM). It was found that the 
effects of grain slots on the downstream flow are 
minimal. Effects dissipate quickly because of rapid 
mixing with the core flow. Carrier et al. [13] used 
an approximate and simple model for particle 
trajectories in a long-bore solid rocket motor. 
Heister and Landsbaum [14] described anomalies 
caused by mass ejection through the nozzle throat of 
Titan 34-D solid propellant rocket motor. Hess et 
al. [15] found that the accumulation of slag in the aft 
end for a short motor is less than accumulation of 
slag for the long motor. Numerical solution of a 
simple-Eulerian-potential - flow / Lagrangian - 
particle - tracking formulation was obtained. It was 
found that gravity acceleration would be discernible 
only for larger particle sizes, above 200-pm 
diameter. Smith-Kent et al. [16] describe a potential 
flow based slag accumulation model with 
Lagrangian particle tracking. Key input parameters, 
namely, particle size, vortex definition, and capture 
criteria are based on empirical data. Johnston et al. 
[17] conducted quasisteady axisymmetric, in viscid, 
coupled two-phase (combustion gas and molten 


AI 2 O 3 ) flow calculations. Six bum back geometries 
(0, 30, 55, 80, 1 10, and 125 sec) that span the total 
bum time of 135 sec, and four different droplet sizes 
(10, 35, 60, and 100 jam) were carried out. They 
have simplified the bimodal distribution droplet size 
of Salita [6] to a flow with just two droplet sizes, 
namely 1.5 and 100 pm. Loh and Chwalowski [18] 
used particles of diameters of 1 to 100 pm in 
converging-diverging nozzles and a mass loading of 
28.8%. Acceleration between lg and 3g has a 
minimal effect on the particles’ behavior in the 
nozzle. Whitesides et al. [19] conducted combined 
analytical and experimental studies to develop an 
understanding of the effects of slag ejection on 
motor performance. A simplistic quasi-steady 
analytical model was formulated for the purpose of 
determining the instantaneous slag flow rate and the 
total quantity of slag required to produce a given 
pressure perturbation. Salita [20] reported slag 
measurements of 58 kg (128 Ibm), 1980 kg (4366 
lbm), and 1102 to 4365 kg (500 to 3500 lbm) in 
SICBM, SRMU, and RSRM, respectively. A flow 
model was proposed for use in predicting slag 
accumulation in these motors. Chauvot et al. [21] 
developed a model to predict slag weight deposited 
in the submerged nozzle of solid rocket motors. In- 
flight acceleration increases slag weight. Sabnis et 
al. [22]conducted two-phase three-dimensional flow 
field in the Titan IV SRM at 17.5 in. bum-back 
geometry using multi-phase Navier-Stokes analysis 
using CELMINT (Combined Eulerian Lagrangian 
Multidimensional Implicit Navier-Stokes Time- 
dependent) code [23]. The geometry in the aft 



closure region is non-axisymmetric due to the nozzle 
being canted with respect to the chamber axis. The 
purpose of the analyses was to understand 
asymmetric insulation erosion in the aft closure 
observed in static and flight tests. Asymmetric 
geometry in the aft closure results in secondary 
flows that can significantly affect the impingement 
pattern of the aluminum oxide droplets on the aft 
closure insulation. Several two-dimensional 
axisymmetric calculations were performed before 
initiating the three-dimensional calculations. They 
were performed to (1) assess the grid resolution, (2) 
obtain consistent flow conditions at the inlet to the 
sixth segment so that the three-dimensional 
simulations could be limited to the region down- 
stream of the fifth segment. The two-dimensional 
axisymmetric analyses comprised segments 1 
through 7, the aft closure propellant and the nozzle 
and solved on two grids. The first grid consisted of 
900 and 90 cells in the axial and radial directions, 
respectively, yielding a total number of 81,000 cells. 
The second grid consisted of 1,175 and 90 cells in 
the axial and radial directions, respectively, yielding 
a total number of 105,750 cells. The viscous sub- 
layer was resolved (i.e., the first point off the wall 
corresponds to y + < 1). Simulation 1 of the present 
study utilized many more cells (Table 1) than in 
Sabnis et al. [22] to obtain y + around 5. This author 
used up to 300 cells in the radial direction in the 
converging-diverging section of the nozzle. The 
size of the cell adjacent to the wall was 10 4 m 
(3.94x1 0’ 3 in.). A smaller size than this would 
generate negative areas which CFD can not handle. 


Therefore, the iterative process between the grid and 
the calculated wall y + ceased. In the three- 
dimensional analyses [22], a grid was used and 
consisted of 351, 90, and 19 in the axial, radial, and 
circumferential directions, respectively, yielding a 
total of 598,500 cells. A propellant combustion 
study conducted by Perkins et al. [24] has shown 
that propellant variability due to subtle raw 
ingredient can affect the quantity and distribution of 
slag formed. Quench bomb tests conducted by 
Brennan [25] at 3.44 MPa (500 psia) chamber 
pressure and particles were quenched at 0.0127 m 
(0.5 in.) from the burning surface. The resulting 
particle size distribution was bi-modal with 50-70% 
(Fig. 10 of Brennan [25] of the particles by weight 
under 5 |im and designated as smoke). The rest of 
the fraction in the 5-700 jam range was designated as 
the coarse fraction. This coarse fraction comprised 
the discrete fraction used in the two-phase CFD 
analysis conducted by Whitesides et al. [26, 27]. 
Whitesides et al. [27] conducted a two-dimensional 
axisymmetric two-phase flow analysis using 
CELMINT code [23]. The overall objective was to 
determine the structure of the flow field in the 
recirculation region underneath the submerged 
nozzle nose and to define the gas flow and particle 
impingement environments along the surface of the 
aft case dome insulation. It was concluded that 
particles were impacting the area underneath the 
nozzle nose and forming a sheet of molten 
aluminum oxide or slag. The sheet flows 
afterwards, along the underneath nozzle nose 
surface as is the direction of the near surface 



velocity vector during the last half of motor bum. 
This slag layer is then sheared from the nozzle 
cowl/boot ring surface and impacts the aft dome 
case insulation at the location of severe erosion. 

Sauvageau et al. [28] discussed the technical 
feasibility of the five-segment reusable solid rocket 
motor approach and quantified the subsequent STS 
capability improvements. The two-dimensional 
axisymmetric analysis (Simulation 1) has the same 
geometry and grid used by Laubacher [29]. The 
analysis in [29] considered adiabatic walls and the 
standard k-e turbulence model solved by a coupled 
solver using in-house code, SHARP* . 

The previous citations are related to internal flow 
in solid rocket motors. They are focused on 

• Internal flow characteristics such as injection 
Reynolds number, turbulence modeling, 
velocity profile, etc. 

• Slag deposition, formation, collection and 
ejection and their relation to chamber pressure 
perturbation and case insulation, and 

• Bore pressure drop. 

• Circumferential pressure gradient in the aft- 
dome of a vectored submerged nozzle. 

• Effects of idealized asymmetric inhibitor stubs 
on circumferential flow. 

• None of the previous studies have dealt with 
convective heat transfer coefficients. 

When considering heat transfer in solid rocket 
motors, surface temperatures and heat fluxes are 
high and very difficult to measure. Ablative 
materials are used to dissipate and inhibit heat 
transfer by erosion and transpiration. It is usually 


estimated using three well-known methods. They 
are the modified Reynolds analogy for laminar flow 
over a flat plate, Dittus-Boelter correlation for fully- 
developed turbulent pipe flow, and the Bartz 
correlation for nozzle flows. 

Related heat transfer studies are given in Refs. [30- 
38]. Bartz [30] extended the well-known Dittus- 
Boelter correlation for turbulent pipe flow to 
account for mass flux and variations in velocity and 
temperature. Back et al. [31 and 33-35] conducted 
analytical and experimental convective heat transfer 
studies in the Jet Propulsion Laboratory (JPL) 
nozzle. Moretti and Kays [32] conducted 
experimental convective heat transfer to an 
essentially constant property turbulent boundary 
layer for various rates of free-stream acceleration. 
Back et al. and Moretti and Kays found that 
acceleration causes a depression in heat transfer rate 
below what would be predicted assuming a 
boundary-layer structure such as obtained for 
constant free-stream velocity. They attributed it to 
re-laminarization of the turbulent boundary layer. 
Moretti and Kays [32] state the above acceleration 
parameter was a result of experimental tests 
conducted in a two-dimensional channel. They 
further state, it is by no means obvious that the same 
acceleration parameter applicable to an 
axisymmetric flow. Wang [37] focused on the 
capability of general-purpose CFD codes in 
predicting convective heat transfer coefficients 
between a fluid and a solid surface. Effects of 
various parameters such as grid resolution, near- 
wall treatments, turbulence models, as well as 



numerical schemes on the accuracy of predicted 
convective heat transfer were studied. Test 
cases included flat plate, pipe flow, JPL nozzle, 
and impinging jets. 

For the lack of reliable thermal conditions, the 
nozzle wall usually assumed to be adiabatic in CFD 
calculations. On the other hand, CFD calculations 
(velocity, density, pressure, temperature, viscosity, 
etc.) and geometry enable someone to calculate heat 
transfer. 

Analyses and Results: 

Two sets of comparative CFD simulations of the 
RSRM and the five-segment concept motor were 
conducted with CFD commercial code FLUENT 
[39, 40]. The first set (Simulation 1) involved a 
two-dimensional axisymmetric model of a full motor 
initial grain RSRM. This analysis was performed to 

• Maintain continuity with previous analyses 

• Serve as a non-vectored baseline 

• Provide a relatively simple checkout tool for 
various CFD solution schemes, grid sensitivity 
studies, turbulence modeling/heat transfer, etc. 
Complete three-dimensional analyses are time 
and cost prohibitive. 

The second set of analyses (Simulation 2) 
included three-dimensional models of the RSRM 
and five segment aft motors with four-degree 
vectored nozzles. The model is composed of a 180- 
degree segment starting at upstream of the aft dome 
and ending at the nozzle exit. This is done to halve 
the number of cells needed. Results of the two 
simulations are detailed in [38] using the turbulence 


Model 1 (Table 1). 

The attributes of the two-dimensional (2D) 
axisymmetric analyses (Simulation 1) are: 

• Significant effort was made to assess grid 

sensitivity and grid consistency with 
turbulence models. 

• Verifying flow/thermal solution quality 

represented by u + vs. y + and t + vs. y + against 
the velocity and thermal laws of the wall in 
turbulent flows. 

• Calculations of nozzle heat transfer, including 

the assessment of turbulent boundary layer re- 
laminarization. Additionally, the viscous 
sublayer thickness is calculated and discussed. 

• Calculations of slag concentration and rate of 

accretion flux and their locations. 

• Calculations of vacuum thrust and specific 

impulse for comparison with the three- 
dimensional (3D) analyses in Simulation 2. 

The attributes of the 3D vectored nozzles 
(Simulation 2) are: 

• Calculations of nozzle heat transfer in vectored 

nozzles. 

• Calculations and presentations of local static 

pressure gradient in the submerged region. 

• Calculations of slag concentration and rate of 

accretion flux sand their locations. 

• Determination and evaluation of 3D effects due 

to nozzle vectoring on vacuum thrust and 
specific impulse in comparison with 2D- 
axisymmetric analysis (Simulation 1). 
Simulation 1 will be discussed first followed by 
Simulation 2. 



Simulation 1: 2D-Axisymmetric 
Analysis for the RSRM at Ignition 

Simulation 1 involves the full motor with initial 
grain configuration. Governing equations, solution 
sensitivity to grid density, boundary conditions, 
computational schemes, numerical convergence 
(residuals), flow fields, turbulence modeling, 
convective heat transfer, and two-phase modeling 
are discussed next. 

Governing Equations: The numerical studies 

considered the solution of the Navier-Stokes 
equations, energy equation, the turbulence kinetic 
energy with its rate of dissipation equations, and the 
necessary constitutive equations (ideal gas law, 
Sutherland’s relation, etc.). The general governing 
equation was 

v.(pv<t> -r, v<p)= s, (D 

and the mass conservation equation 

V. (p v) = 0 < 2 > 

where <f> can be velocity components («, v, n ), 
enthalpy (/), turbulence quantities (k, e), and species 
concentrations; r is an exchange coefficient for 0; 
S 0 is a source term for 0 per unit volume. 

Solution Sensitivity to Grid Density: Grids used 
in Simulations 1 and 2 are summarized in Table 1. 
The grids used in the two-dimensional axisymmetric 
analyses (Simulation 1) are quadrilateral or 
triangular. The grids in the two-dimensional study 
were designed, solved, and iterated on to give the 
desired values for y + so that consistency with the 
turbulence models was achieved as given in Table 1 . 


The grids used in the three-dimensional are 
hexahedral and/or tetrahedral. Grid adaption was 
used in some grids as given in Table 1. Adaption 
worked well in terms of generating similar flow 
fields (Mach number and pressure distributions) and 
lowering the wall v + to the desired value. On the 
other hand, the calculated wall y + and h/Cp profiles 
have no distinct trend or profile as known from 
measurements and approximate methods in nozzle 
flows. Therefore, it was discontinued in Grid 4 of 
Table 1. All the grids were generated by using 
GRIDGEN [41] and made orthogonal and smoothed 
from one domain to another in the case of quads and 
hexahedral. 

For simple geometries , quadrilateral or 
hexahedral meshes can provide high quality 
solutions with fewer cells than comparable 
triangular or tetrahedral meshes [39]. For complex 
geometries , quadrilateral/hexahedral meshes show 
no numerical advantage, and time can be saved 
using triangular/tetrahedral mesh. This would 
enable refined grids in regions of interest along with 
the consistent turbulence model using adaption 
techniques based on wall v + . The convective heat 
transfer coefficient and wall v + are analogous and 
interdependent as will be shown. 

Boundary Conditions: the boundary conditions in 
Simulation 1 are as follows: 

At the propellant surface : mass flux is calculated 
as a function of the local static pressure as 

m =«, P p a [P (*, y)] n 

In addition, uniform chamber temperature, flow 
direction that is normal to the propellant surface, an 
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assumed turbulence intensity [36], (/), and hydraulic 
diameter (dh) were specified. 

References 42, 43 and 44 give the chamber 
pressure as 6.30 MPa (913.85 psia), 6.24 MPa (905 
psia), and 6.27 MPa (910 psia), respectively. 

The total pressure used in this study is 6.25 MPa 
(906 psia). This chamber pressure along with 
propellant formulation was used as input to the 
NASA-Lewis program [45] or the ODE module of 
the SPP code [42] to obtain chamber temperature, 
gas molecular weight, dynamic viscosity, thermal 
conductivity, and specific heat at constant pressure. 
They are 3416.04 K (6148.87 R), 28.482 kg/kgmole, 
9.643x10 5 N-s/m 2 (6.48x1 0' 5 Ibm/ft-sec), 0.383 
W/m-K (0.221 Btu/hr-ft-R) and 2005.48 J/kg-K 
(0.479 Btu/lbm-R), respectively. The above 
properties yield a value for Pr of 0.51. Propellant 
density is taken as 1743.83 kg/m 3 (0.063 lbm/in 3 ). 
The reference pressure is taken as 4.30 MPa (625 
psia). The reference bum rate and exponent are 
taken as 9.347x1 0' 3 m/s (0.368 in/s) and 0.35, 
respectively. The above parameters yield a value 
for a in Eq. (3) of 0.07765 kg 065 m 165 s "° \ 
Turbulence intensity of 5% [36] was specified. 

At wall: No slip condition was specified as 
velocity boundary condition. Isothermal wall was 
specified (3,000K (5,400R)) at the submerged wall. 
This is because calculations show the temperature to 
be uniform and the velocity is small in comparison 
with the accelerating flow in the nozzle. On the 
other hand, a specified axial surface temperature 
was assumed for the converging-diverging part of 
the nozzle and will be given later in the heat transfer 


section. 

At exit: Supersonic boundary condition where 
the quantities (P, 7, w, v, k, e) are calculated from 
cells upstream of the exit was specified. 
Computational Schemes: The segregated solver in 
the commercial code Fluent [40] was used. 
Differencing schemes utilized are 1 st and 2 nd order 
Upwind, Power law, and Quick schemes. The 1 st 
Upwind scheme was used to start the problem and 
the higher order schemes to obtain the final results. 
The 2 nd order Upwind and Quick schemes were 
found to give similar results in terms of mass flow 
rate and mass imbalance, head-end pressure, 
chamber pressure drop, and maximum Mach number 
at the nozzle exit. 

Numerical Convergence (Residuals): Numerical 

convergence was achieved by satisfying four 
requirements in the following sequence. First, the 
residual error diminished as the number of iterations 
was increased. Second, the profiles of the variables 
ceased to change, at least qualitatively. Third, a first 
monitor on the total pressure at the propellant 
surface or at the flow inlet until the average total 
pressure ceased to change. Fourth, a second monitor 
on the mass imbalance between the inlet (propellant 
surface) and the outlet (nozzle exit) mass flow rates 
until it reached a small value ( 10 3 - 10' 5 kg/s). 
Flow/Thermal Fields: Figure la shows the local 
Mach number distribution in the full motor. The 
maximum Mach number obtained is 3.56 at the 
nozzle exit plane at the centerline. Similarly, Figs. 
1 b( 1 -5) shows the calculated Mach number in the 
cavity of the submerged region using the four 



turbulence Models (Table 1 ). They are plotted for 
the same range (0 < M > 0.01) and the same scale 
(30,000) except Model 1 with the undersigned grid 
for y + where the scale was reduced to a 1 ,000. As 
the model is refined with the appropriate turbulence 
model, the re-circulation profile changed. Model 1 
with the undersigned grid over-predicts the 
circulation motion. Models 1 and 2 are similar. 
Similarly, Models 3 and 4 are similar. Obviously, 
these models would have different effect on slag 
concentration and rate of accretion flux and will be 
discussed in relation to the discrete phase 
concentration and rate of accretion flux. 

Figure lc show the static temperature and stream 
function, respectively. A drop of 1750K (3,150 R) 
in static temperature has occurred. 

Bore Pressure Drop: Figure 2a shows the static 
pressure distribution in the whole motor. Figure 2b 
shows the local axial static pressure along the 
centerline of the RSRM. The interest here is to 
compare the calculated pressure drop against 
measurements, matching head-end total pressure and 
numerical accuracy. Four solutions were conducted 
using a segregated solver along with a 2 nd order- 
Upwind differencing and in conjunction with the 
four turbulence models given in Table 2. A fifth 
solution was conducted using a coupled solver with 
the first turbulence model. The viscosity ratio 
(effective (laminar plus turbulent) to laminar) was 
limited to 10,000 to match the total pressure of 6.25 
MPa (906 psia). Without limiting this viscosity 
ratio to 10,000, the calculated chamber pressure (at 
some locations) is unrealistically higher than the 


known total chamber pressure with some 
differencing schemes. 

Also shown the measurements of the local static 
pressure along QM-7 and QM-8 test motors. 
Qualitative agreement of the calculated local static 
pressure is achieved. 

The ratio of the calculated head end pressure to 
the specified head end pressure (P (K max / P (t ) is 
calculated as 1.02, 1.02, 1.02, 1.02, and 1.06, 
respectively. Clearly, the head end pressure is over- 
predicted by 2 percent and 6 percent using the 
segregated and the coupled solver, respectively. 
Upstream of the nozzle throat, the ratio of the 
calculated head end pressure to the specified head 
end pressure (P t >(x/r* = 64, y = 0) / P (> ) is calculated 
as 0.8, 0.80, 0.80, 0.80, and 0.87, respectively. 
Clearly, the head end pressure is over-predicted by 2 
percent and 6 percent using the segregated and the 
coupled solver, respectively. Thus, the 
corresponding non-dimensional chamber pressure 
drops are 0.22, 0.22, 0.22, 0.22, and 0.19, 
respectively. The first four solutions match the 
head-end pressure but not the pressure drop. The 
fifth solution over-predicts the head-end pressure 
but matches the chamber pressure drop. The 
recommended pressure drop is 18 percent as given 
from [29, 44]. It is to be noted that the solution in 
[29] was conducted on Grid 1 (68,713 cells) with 
adiabatic walls, coupled solver (SHARP code), 
turbulence intensity of 5% and viscosity ratio of 
10,000. 

Turbulence Modeling: Four turbulence models are 



(4b) 


utilized in this study and given in Table 2. Figure 3a 
shows the calculated wall y + along the converging- 
diverging part of the nozzle using the four 
turbulence models. The non-dimensional velocity in 
wall coordinates (w + ), non-dimensional distance 
from wall in wall coordinates (y + ) and frictional 
velocity are defined [40, 46] as 

+ u(y) + “ r (y w -y) (4a) 
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The first two models give similar profiles and values 
from Grid 5. Since the values are between 20 and 
100, the first and second turbulence models used are 
consistent with grid resolution. The third and fourth 
turbulence models give similar profiles and values 
from Grid 6. The wall y + calculated values using 
Grid 6 are shown to be below 8; therefore, using the 
third and fourth turbulence models are consistent 
with grid resolution. The results of calculated wall 
y + for the third and fourth turbulence models are 
also shown separately in Fig. 3b for clarity. The 
calculated v + could not be dropped without using a 
lot more cells in the stream-wise direction. The 
minimum spacing was specified as small as 10' 4 m 
(3.94x10’ 4 in.). 

Figure 3c shows the velocity law of the wall 
behavior of the boundary layer at the plane of the 
nozzle exit. The results using four turbulence 
models are shown. Also shown is the Spalding’s 
velocity law of the wall [46] for incompressible 
flow: 
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where k and B are taken [46] as 0.4 and 5.5, 
respectively. 

The first two models fall in the log-law and wake 
regions. The first flow cell is located at y + of 22. 
The third and fourth turbulence models compare 
well with the law of the wall in the inner layer but 
not in log-law region. The first flow cell is located 
at y + of 1.5. Similarly, Fig. 3d shows the 
temperature law of the wall at the plane of the 
nozzle exit. The temperature law of the wall is 
given [47] as 


t + = Pr y + ;y + <13.2 
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where the non-dimensional temperature in wall 
coordinates is defined [47] as 
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and where Pr, and / are taken as 0.85 and 0.41, 
respectively. The results using four turbulence 
models are shown. They fall on top of each other 
but they do not extend the same values of v + . Also 
shown is the thermal law of the wall for air taken 
from Kays and Crawford [47] for turbulent flow. 
Convective Heat Transfer: Figure 4a shows the 
submerged and nozzle walls and the specified 



surface temperature used as thermal boundary 
conditions. The submerged wall was assumed to be 
isothermal at 3,OOOK (5,500R). On the other hand, a 
specified axial surface temperature was assumed for 
the converging-diverging wall of the nozzle as 
shown in Fig. 4a. The surface temperature profile 
was taken from [48] and curve-fitted by this author 
using TableCirve2D [49] as 

_ a +cx+ex 2 + g x* +/ x 4 +kx* O) 

1 + b x + d x 2 + / jc 3 + h x 4 + j x 5 

where x is taken along the nozzle surface and where 
the coefficients are given as follows: a = 2789.03, b 
= - 4.61, c = -13307.21, d = 1 1.96, e = 32905.45,/ = 
- 6.33, g = - 20712.33, h = - 0.916,; = 1 1 74.06 and 
k = 815.51. The correlation coefficient, r 2 = 0.997. 
A user defined function (UDF) was used to compile 
the specified surface temperature profile. Assuming 
this surface temperature profile enables the 
calculation of the heat flux that in turn enables the 
calculation of the convective heat transfer 
coefficient depending on the assumption of a 
reference temperature. Calculated centerline and 
recovery temperatures are also shown and will be 
discussed shortly. 

In the present study, the ratio of the convective 
heat transfer coefficient and the specific heat at 
constant pressure is calculated by two methods. 
This ratio is used as an input in the CMA code [50] 
for nozzle erosion predictions. In Method 1, it is 
calculated internally using the calculated heat flux 
based on the difference between the local specified 
surface temperature and the chamber temperature 
(used as a reference temperature for simplicity), i.e.. 
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and shown in Fig. 4b(l). Figure 4b(2) is a 
magnification of the throat region and taken from 
Fig. 4b(l) for clarity 

In Method 2, it is calculated using the recovery 
temperature and defined as 

(7a) 
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And shown in Fig. 4b(3). The recovery temperature 
[47] is given as 

= r (x)+^[7;, -Tjx)] (7b) 
and where 9t = Pr m is used for turbulent flow and is 
calculated to be 0.8. The temperature at the edge of 
the boundary conditions has been replaced by the 
local axial static temperature along the motor 
centerline and shown Fig. 4a. The recovery 
temperature is calculated at the motor centerline and 
shown in Fig. 4a to be less than the chamber 
pressure. The chamber and recovery temperature 
are matched only by taking Pr as 1 which yields I 
for 9? which correspond to an ideal situation. Figure 
4b(4) is a magnification of the throat region and 
taken from Fig. 4b(3) for clarity. 

The above two heat transfer coefficients are 
calculated in conjunction with the four turbulence 
models given in Table 2. Also shown, is the 
calculated h/C P from measured data [48] to match 
maximum measured erosion. This measured h/Cp is 
a result of heat transfer (convective, radiative and 
conductive), chemical, and particle impingement. 

Furthermore, additional results from the 
Turbulent Boundary Layer (TBL) code [51] are also 



shown. The TBL predictions overestimate the ratio 
of h/Cp in most part of the converging-diverging 
part of the nozzle. The TBL routine [51] of the SPP 
code [42] is an integral method that uses the 
boundary-layer edge conditions from the two- 
dimensional, two-phase flow program of SPP 
program [42]. On the other hand, it underestimates 
the ratio in the last one-third of the diverging section 
of the nozzle. This may be due to the lack of 
particle impingement treatment in the TBL code 

[51]. 

The measured values [48] of h/Cp are the net 
result of many effects such as convective, radiative, 
particle impingement, chemical reactions, and 
others. Therefore, calculating h/Cp from Methods l 
and 2 and comparing them with the measured data 
[48] in a relative manner would be a sound 
approach. The following is a discussion of the 
results of each turbulence model used. 

1A) Based on Method 1 ( Eq. (6Lk 

(1) Model /, Grids 2 and 5: The first turbulence 
model on Grids 2 and 5 overestimates the ratio of 
h/Cp in the whole nozzle except in a small portion at 
the nozzle exit. It underestimates it in the diverging 
section of the nozzle. The ratio h/Cp drops at the 
nozzle throat and then increases then decreases in 
the convergent section of the nozzle. 

(2) Model 2, Grid 5: The second turbulence 

model matches the measured ratio of h/Cp in the 
converging part of the nozzle. It drops significantly 
at the throat and increases again and underestimates 
it in the diverging section of the nozzle. Again, the 
calculated ratio h/Cp does not include the effect of 


particle impingement. 

(3) Model 3> Grid 6: The third turbulence model 
matches the measured ratio of h/Cp in the 
converging part of the nozzle. It underestimates it at 
the throat region. It underestimates it in the 
diverging section of the nozzle. Again, the 
calculated ratio h/Cp does not include the effect of 
particle impingement. 

(4) Model 4 , Grid 6: The fourth turbulence 

model behaves similarly to the third model. 

It can be concluded that Model 2 on Grid 5 
gives the best results in relative comparison with the 
measured data [48] as shown in Fig. 4b(2). 
Furthermore, it does not require the grid refinement 
as required for the two-layer zonal model. Thus the 
execution time is much less. 

IB) Based on Method 2 (Eq. (7)): 

The results of Method 2 differ from Method 1 in 
the divergent section of the nozzle. Method 2 is 
more sensitive at the throat and downstream at the 
throat than Method 1. It overestimates the ratio 
h/Cp at the throat. Method 2 neither follows the 
profile qualitatively from the measured values nor 
follows the usual profile from Bartz correlation [30]. 
The author concludes the Method 1 is more 
appropriated than Method 2. 

The four turbulence models show three common 
behaviors: 

• First: the maximum convective heat transfer 
coefficient is upstream of the nozzle throat. 

• Second: the profile of the calculated ratio h/Cp 
from Simulation 1 and along the nozzle wall is 
analogous to the y + distribution shown in Fig. 3a 



and 3b. This confirms the interdependency of 
heat transfer on the turbulence model used, is 
directly consistent with the grid resolution. 

• Third: in the vicinity downstream of the nozzle 
throat, the ratio h/Cp decreases suddenly and 
increases again. 

The sudden drop in h/Cp is attributed to the large 
drop in the specified surface temperature (Fig. 4a). 
In the vicinity of the throat, the surface temperature 
dropped by 590K within 0.5 m (1062 R within 1.64 
ft). This was easily verified by specifying uniform 
surface temperature at 2000 K, respectively. The 
sudden drop in h/Cp was not calculated using an 
isothermal boundary condition. The nozzle wall is 
an ablating (transpiring) surface. The boundary 
layer adjacent to the nozzle wall would be oxidizer- 
poor that inhibit chemical reactions and heat 
transfer. A transpiring surface is not modeled in this 
study. It may show a drop in the convective heat 
transfer. 

Figure 4c shows the local shear stress along the 
nozzle wall calculated using four turbulence models 
as given in Table 2. Figure 4d shows the 

corresponding local skin friction coefficient along 
the nozzle wall. It has the following usual 
definition. 

C f _ T s (8) 

2 

Figure 4e shows the corresponding local Reynolds 
number along the centerline based on local nozzle 
diameter. 

For fully developed pipe flow, the non- 
dimensional pressure gradient along the nozzle wall 


is calculated as 
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and shown in Fig. 4f. The Van Driest constant is a 
measure of viscous sublayer thickness and is given 
[47] as 

26 (9b) 
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where a , b, and c are taken as 7.1, 4.25, and 10.0 
[47], respectively. It is equal to 26 for a zero 
pressure gradient boundary layer without wall 
transpiration. Assuming no wall transpiration, A + 
was calculated and shown in Fig. 4g. It has the 
opposite profile of P/ + . 

The corresponding Pi + and A:+ for a simple 
external boundary layer are calculated respectively 
as 

K (9c) 
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and shown in Figs. 4h and 4i. In the above, K is an 
acceleration parameter and is calculated from [47] 

as 
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and shown in Fig. 4j. The calculated values of K are 
smaller then the transition value of 3xl0 6 . 
Therefore, re-laminarization of the turbulent 
boundary layer did not occur. Additionally, the non- 



dimensional pressure gradient is defined as 

fl(dP/dx) <90 
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Two-Phase Modeling: The two phases are 
continuous phase (combustion gas) and a discrete 
phase (aluminum oxide droplets (AI 2 O 3 )). Since 
droplet agglomeration/breakup is not modeled, they 
are assumed spherical and of a constant size 
distribution anywhere in the motor. Solid spherical 
particles are injected at propellant surface 
(Simulation 1) and at the flow inlet (Simulation 2). 

The particle size and distribution used in this 
study were based on investigation [25] into the 
effects of ingredient property variation on RSRM 
TP-H1148 production propellant. Aluminum oxide 
caps fraction is the total fraction of aluminum that is 
converted to large Al 2 Oi discrete phase particles 
after combustion is complete and is taken from [27] 
as: 


f AI = 0.2833 ( 10a) 

J Alu min a 

The remaining aluminum oxide in the form of 
smoke is therefore 


f =1 - f . =0.7167 (10b) 

J smoke J A/u minu 

The discrete phase portion of the total particulate 
population (Alumina) is taken as 

f DP = f Alu = 0.0944 ( ' ° C ' 

The mass flow rate of the discrete phase is 
calculated in relation with the continuous phase 
(gas) as follows: 
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Inert particles of density of 1601.85 kg/m 3 (100 
lbm/ft 3 ), thermal conductivity of 10.47 W/m-K (6.05 
Btu/hr-ft-R), and specific heat of 1800.32 J/kg-K 
(0.43 Btu/lbm-R) are specified. The molecular 
weight of the gas was calculated as 

18.76 + 7.58 / . (10e) 

MW - J smoke 

0.926 

which gives 26.07 kg/kgmole. 

Particle size distribution for the alumina included 
in the discrete phase was defined by fitting the size 
distribution to the Rosin-Rammler equation [Chap. 
14 of 40]. In this approach, the complete range of 
particle sizes is divided into a set of discrete size 
ranges, each to be defined by a single stream that is 
part of the group. The particle size data obeys the 
distribution shown in Fig. 5a and is curve-fitted as 




a exp 



(100 


where the coefficients a and b are calculated to be 
1.01475 and 86.5651, respectively and with a 
coefficient of determination of 0.99. 

The Rosin-Rammler distribution function is 
based on the assumption that an exponential 
relationship exists between the particle diameter, 
(D), and the mass fraction of particles with diameter 
greater than D, Mo as 


M D = exp 
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where D is the mean diameter and N is the spread 
parameter. The average diameter is obtained by 
noting that the value of D at which Mo = exp(-l) = 
0.368. From the above figure, the average diameter 
is 87.53 pm. The numerical value of the spread 
parameter is calculated from the above relationship 



as 


N In (- In M 0 ) (HO 

Substituting the values, an average value for N was 
calculated to be 1.1426. In summary, 8 discrete 
particle sizes ranging from 10 to 600 microns with a 
mean diameter of 87.80 jim and a spread of 1.1428 
were specified. 

In Simulation 1, particles are injected at zero 
velocity and at 0% temperature lag (equal to gas 
temperature). The walls of the submerged region 
and nozzle are set as trap walls for particles upon 
impact as boundary conditions. Propellant surface 
was set as a reflecting surface for particles. 

Particle concentration is a flow quantity that 
exists in the flow field and at the wall. The particle 
trajectories and their local A1 2 Oj concentrations are 
of interest to understand the submerged region and 
nozzle environments in terms of insulation 
performance and possibly Al 2 Ch slag rate of 
accretion flux or accumulation. Figures 5b and 5c 
show the calculated Al 2 Ch local concentration along 
the submerged and nozzle walls, respectively. The 
maximum calculated concentrations are found to be 
6.28 kg/nv (0.39 lbm/ft 3 ) and 0.3 kg/rrv (0.02 
lbm/fr ), respectively. 

Figures 5d and 5e show the calculated Al 2 Ot rates 
of accretion fluxes along the submerged and nozzle 
wall, respectively. The maximum calculated rate of 
accretion fluxes are found to be 15 kg/m 2 -s (3.07 
lbm/fr-s) and 12 kg/rrf-s (2.46 lbm/ft 2 -s), 
respectively. The maximum occurs further upstream 
on the submerged wall. 


Simulation 2: Four-degree Vectored 
Nozzles for the RSRM and FSM 

Results similar to those for Simulation I are 
sought in Simulation 2. Additionally, the effects of 
nozzle vectoring on the above parameters are also 
included and discussed. 

Simulation 2 involves two three-dimensional aft 
motors with vectored nozzles for the RSRM and the 
FSM as shown Figs. 6a and 6b, respectively. The 
solution domain extends from 0 to 180 degrees 
segment with a plane of symmetry boundary 
conditions were used for the 0 to 180 degrees. 
Table 1 gives motor chamber, nozzle throat, and 
nozzle exit diameters and the corresponding area 
ratios. Differences in motor geometries are not clear 
in these two figures. Therefore, Figs. 6c and 6d are 
shown to give the details of the 3D grid on the 
symmetry planes of the RSRM and the FSM motors, 
respectively. The two nozzles are vectored 
downward. This causes the bottom halves of the 
RSRM and FSM nozzles, to protrude ahead in the 
flow field. In addition, the cavity at the bottom 
becomes larger than the cavity on the top side. 
Similarity between the grids will be verified by the 
calculated wall y + and will be given later. 

Grid-independent results in the three-dimensional 
analyses are cost and time prohibitive. In addition, 
since the interest is in comparing the results of the 
two motors (in terms of an increase or decrease), the 
approach was taken to make the two grids as 
comparable as possible in terms of similar number 
of cells (approximately 325,000 hexahedral cells 
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each) and similar local clustering of cell sizes in the 
gridding process. A finer grid for the FSM motor 
was conducted for sensitivity studies and given in 
Table 1. 

Boundary Conditions: The geometry considered is 
a representative of late bum time. The boundary 
conditions in Simulation 2 are as follows: 

At the inlet: An assumed uniform total pressure 
and temperature, calculated velocity that is normal 
to inlet, an assumed turbulence intensity (/), and a 
hydraulic diameter (dh) were specified 

P n =c,T a =c,e = o on 
A total pressure (4.25 MPa (616 psia)) and total 
temperature (3368 K (6062.4 R)), turbulent intensity 
of 5% and the hydraulic diameters (Table 1) were 
specified. The velocity components are calculated. 
This chamber pressure along with propellant 
formulation was used as input to the NASA-Lewis 
program [45] or the ODE module of the SPP code 
[42] to obtain chamber gas temperature, molecular 
weight, dynamic viscosity, thermal conductivity, and 
specific heat at constant pressure. They are 
respectively 3368.32 K (6062.96 R), 28.378 

kg/kgmole, 9.582x10° N-s/m 2 (6.439x10 5 lbm/ft- 
sec), 0.383 W/m-K (0.22 Btu/hr-ft-R) and 2002.97 
J/kg-K (0.478 Btu/lbm-R). The reference pressure 
is taken as 4.3 MPa (625 psia). Turbulent intensity 
of 5% [36] and the hydraulic diameter were 
specified. 

At wall: No slip condition for the velocity. The 
submerged wall was assumed to be isothermal at 
3,000K (5,500R). On the other hand, a specified 
axial surface temperature was assumed for the 


converging-diverging wall of the nozzle as shown in 
Fig. 4a as in Simulation 1 . 

At exit: Supersonic boundary condition where 
the quantities (P, 7, w, v, w\ k, e) are calculated from 
cells upstream of the exit. 

Computational Schemes: A segregated solver was 
used. Differencing schemes utilized are 1 st and 2 nd 
order Upwind, Power law, and Quick schemes. The 
1 st Upwind scheme was used to start the problem 
and the higher order schemes to obtain the final 
results. The 2 nd order Upwind and Quick schemes 
were found to give similar results in terms of mass 
flow rate and mass imbalance, head-end pressure, 
chamber pressure drop, and maximum Mach number 
at the nozzle exit. 

Numerical Convergence (Residuals): Numerical 
convergence was achieved by satisfying four 
requirements in the following sequence. First, the 
residual error diminished as the number of iterations 
was increased. Second, the profiles of the variables 
ceased to change, at least qualitatively. Third, a first 
monitor on the total pressure at the propellant 
surface or at the flow inlet until the average total 
pressure ceased to change. Fourth, a second monitor 
on the mass imbalance between the inlet (propellant 
surface) and the outlet (nozzle exit) mass flow rates 
until it reached a small value (10’ 3 - 10° kg/s). 

Flow Fields: Figures 7a and 7b show the flow field 
depicted in terms of the velocity magnitude in the 
RSRM and the FSM motors respectively. A 
circumferential flow developed from the bottom half 
(due to a circumferential pressure gradient) to the 
upper half (lower pressure) in the submerged 



regions of the two motors. The circumferential flow 
emanates from the bottom half and moves to the 
upper half. This means that, under vectoring, the 
upper and lower halves experience different 
flow/thermal environments. The flow originates at 
bottom, moves circumferentially, and is flushed 
from the top half. 

Figures 7c and 7d show the velocity magnitude in 
the aft motor chamber for the vectored nozzle 
simulations of the RSRM and the FSM, 
respectively. Average inlet velocities are calculated 
to be 1 16 m/s (380.6 ft/s) and 93 m/s (305.1 ft/s) in 
the FSM and RSRM, respectively. The resulting 
ratio between the inlet velocity magnitudes is 1.25. 
In the case where the chamber pressure and 
temperatures are the same as in this analysis, this 
translates into an increase of 25% in the mass flow 
rate of the FSM in comparison with the RSRM. 
This is only due to the opening of the throat in the 
FSM in comparison to the RSRM. Vortices have 
developed in the bottom halves of the two motors 
and not in the top. 

Figures 8a and 8b show the Mach number 
distributions in the vectored nozzles of the RSRM 
and the FSM, respectively. The overall profiles in 
the two motors are qualitatively similar. The 
maximum Mach numbers were calculated to be 3.33 
and 3.25 in the RSRM and FSM, respectively. They 
are respectively 93% and 91% of the maximum 
Mach number (3.54) calculated in the 2-D 
axisymmetric (Simulation 1 ). The FSM geometry 
has a larger throat and a smaller expansion area ratio 
in comparison with the RSRM. Grid 3 of the FSM 


comprises 1,129,104 hexahedral and uniform cells 
which are approximately 3.5 times larger than the 
cells in Grid 1 (324,000 hexahedral cells). It is not 
designed to lower wall v + as in other grids. Rather, 
it is to verify sensitivity of the decrement in 
maximum exit Mach number. The maximum 
calculated Mach number on a larger grid (Grid 3) 
for the FSM is calculated to be 3.25. Therefore, the 
difference in exit Mach numbers between the 2D- 
axisymmetric and 3D appears to be because of 
vectoring and not the grid size or refinement. 

Figures 9a and 9b show the local static pressure 
distributions in the vectored nozzle simulations of 
the RSRM and the FSM, respectively. The overall 
profiles are qualitatively similar. The total 
boundary condition was set at the inlet of both 
motors to be the same. As expected, the overall 
static pressure distribution is the same in the two 
chambers. The minimum static pressure at the exit 
of the RSRM (0.23 atm (3.39 psia)) is lower than 
the minimum static pressure at the exit of the FSM 
(0.31 atm (4.52 psia)). This is because of the slight 
difference in expansion area ratio. The static 
pressure is approximately equal to the total pressure 
set at the flow inlet boundary condition, which is 
approximately 4.25 MPa (616 psia). 

Figures 9c and 9d show the local circumferential 
static pressure distributions in the vectored nozzles 
of the RSRM and the FSM, respectively. The planar 
distributions are selected in the submerged region of 
each motor. They are geometrically parallel to the 
flow inlet and perpendicular to the flow direction. 
They show the pressure to vary axially and 



circumferentially. The variation increases as the 
axial distance increases in the submerged region. 
The max variation of static pressure in the cavity 
volume amounts to 0.14 to 0.27 atm (2 to 4 psid). 
The circumferential static pressure in the deepest 
plane is a fraction of this maximum. 

Pressure Gradient in the Submerged Region of 
the Nozzle: Figures 9c and 9d show the static 
pressure gradient in the submerged regions of the 
nozzles of the RSRM and FSM. It varies axially, 
radially and circumferentially. Figure 9e is a 
simplified version. The deepest plane has been 
reduced to a ring to narrow the circumferential 
pressure gradient. 

In the axial direction : In the RSRM, the present 
results show a maximum axial drop in the static 
pressure (at 0 = 40°) of 0.06 atm (0.87 psid). In the 
FSM, the corresponding maximum axial drop in the 
static pressure (at 0 = 40°) of 0.1 1 atm (1.68 psid) 
was calculated. 

In the radial direction: A maximum of 0.01 atm 
(0.1 psid) and 0.01 atm (0.18 psid) circumferential 
pressure gradient occurred in the RSRM and the 
FSM, respectively. 

In the circumferential direction: In the RSRM, 
the present results show a maximum axial drop in 
the static pressure (between 0 = -90° and 40°) of 0.0 
atm (0.02 psid). In the FSM, the corresponding 
maximum axial drop in the static pressure (between 
0 = -90° and 40°) of 0.07 atm (1.0 psid) was 
calculated. For engineering design, the 
circumferential pressure differential between the 
chamber pressure at the angle of 0 = -90° to the 


angle of 0 = 30° should be considered. This would 
result in a maximum circumferential pressure 
differential of 0.06 atm (0.87 psid) and 0.11 atm 
(1.68 psid) in the RSRM and FSM, respectively. 

Also shown in the figure, cold flow experimental 
results from Whitesides et al [7]. The tests were 
conducted in the NASA Marshall Space Flight 
Center Air flow facility with a 7.5 percent scale 
model of the aft segment of the RSRM. Static and 
dynamic pressure were measured in the model for 
nozzle vector angles of 0°, 3.5°, and 7° with 
propellant configuration corresponding to motor 
bum times of 0, 9, 19, and 114 s. The experimental 
results [7] as given assume no axial pressure drop. 
The maximum circumferential pressure gradient 
occurs at 19 sec and amounts to 0.07 atm (1.05 
psid). This is larger than the calculated value in the 
present calculation. 

Turbulence Modeling: Figures 10a and 10b show 
the local calculated wall y + of the RSRM and the 
FSM, respectively. Two curves are shown for the 
submerged wall and 2 curves for the nozzle wall. 
One curve pertains to the bottom surface and the 
second pertains to the top surface and 180 degrees 
apart. The two surfaces lie on the symmetry plane. 
Note the lower half protrudes slightly ahead of the 
upper half. Since most of the calculated wall v + are 
higher than 100 on this coarse grid, the standard k-£ 
model with the standard wall functions (high Re k-e 
turbulence model) and k-£ RNG (renormalization 
group theory) with standard wall functions were 
used to calculate the final results on Grid 1. The 
similarity of the calculated y + values for the two 
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motors indicates the similarity of the two grids. 
This is not fortuitous or accidental. Two 
simultaneous GRIDGEN sessions were active. 
Adjustments made (number of cells and spacing) for 
one grid were duplicated on the other and vice versa. 
The wall y + is a function of Mach number and cell 
size adjacent to the wall. The maximum calculated 
wall v + in the submerged regions of the RSRM and 
the FSM are 1871 and 1986, respectively. 
Similarly, the maximum calculated wall v + in the 
nozzle regions of the RSRM and the FSM are 1861 
and 1852, respectively. 

Heat Transfer: The specified surface temperatures 
at the submerged and nozzle walls used in the two- 
dimensional axisymmetric analysis (Simulation 1) 
are used here in the 3D vectored nozzle of the 
RSRM and FSM, respectively. It is shown in Fig. 
4a and calculated as a function of only the axial 
location for each motor. It is invariant in the 
circumferential direction. 

The calculated profiles are shown in Figs. 11a 
and 1 lb for the RSRM and FSM, respectively. The 
maximum calculated h/Cp in the submerged regions 
of the RSRM and the FSM are 4.59 kg/rrf-s (0.94 
lbm/fr-s) and 4.54 kg/m 2 -s (0.93 lbm/fr-s), 
respectively. Comparable convective heat transfer 
environment in the chambers of the RSRM and the 
FSM have been calculated. Similarly, the maximum 
calculated Ir/Cp in the nozzle regions of the RSRM 
and the FSM are 4.20 kg/m -s (0.86 Ibm/ft'-s) and 
4.74 kg/m 2 -s (0.97 lbm/fr-s), respectively. An 
increase in the convective heat transfer environment 
of 13% in the FSM in comparison of the RSRM has 


occurred. 

Two-Phase Modeling: Same assumptions and 

treatments used in Simulation 1 were used in 
Simulation 2 with one exception. In Simulation 2, 
particles are injected at the flow inlet (vs. propellant 
surface as done in Simulation 1). In both 
simulations, particles are injected at approximately 
zero velocity lag and zero temperature lag. This is 
because at that axial location (flow inlet), it is 
assumed the particles have accelerated and assumed 
the gas velocity with very small velocity lag. 

Figures 12a and 12b show the local discrete 
phase model (DPM) of AI 2 O 3 concentration 
calculated in the RSRM and FSM, respectively. It 
shows local variation between the top and bottom 
sides in each motor and between the two motors. 
The maximum AI 2 O 3 concentrations in the 
submerged regions of the RSRM and FSM were 
calculated to be 6.51 kg/m 3 (0.41 lbm/ft 3 ) and 3.37 
(0.21 lbm/ft 3 ) in the RSRM and FSM, respectively. 
The maximum AhO< concentration occurred at the 
converging part of nozzle and on the half that is 
protruding ahead in the flow due to vectoring. The 
maximum AI 2 O 3 concentrations in the nozzle regions 
of the RSRM and FSM were calculated to be 0.35 
kg/m' (0.02 lbm/ft 3 ) and 0.43 (0.03 lbm/ft 3 ) in the 
RSRM and FSM, respectively. 

Figures 13a and 1 3b show the discrete phase 
model (DPM) AI 2 O 3 rates of accretion fluxes in the 
RSRM and FSM, respectively. The slag rates of 
accretion fluxes in the submerged regions of the 
RSRM and FSM were calculated to be 41 .52 kg/rrf- 
s (8.50 lbm/fr-s) and 42.2 kg/m“-s (8.64 lbm/fr-s). 



respectively. Similarly, slag rates of accretion 
fluxes in the nozzle regions of the RSRM and FSM 
were calculated to be 6.50 kg/m 2 -s (1.33 Ibm/fr-s) 
and 7.6 kg/m“-s (1.56 lbm/ft -s), respectively. A 
similar rate of accretion flux occurred in the FSM 
and in the RSRM. 

It is to be noted that no particle break-up was 
modeled. Therefore, the calculated concentrations 
or rates of accretion fluxes may not be accurate in 
the nozzle past the throat. They are meant to show 
the trend. 

Vacuum Thrust and Specific Impulse 
Calculations: The nozzle thrust was calculated [52] 
as 

F = jj u{pV .n)dA + {P t -PjA t ( |2a ) 

and vacuum thrust was calculated as 

F vac =$u(pV.r,)dA + P e A e d2b) 

and the vacuum h p is calculated as 

F (12c) 

j vac 

1 Sp , V'£Jt ~~ 

m g 

A user defined function (UDF) was written to 
calculate mass flow rate, thrust, and vacuum I sp on 
demand. Chamber pressures and nozzle expansion 
ratios used are shown in Figs. 14a and 14b, 
respectively. The calculated mass flow rates, 
vacuum thrust and vacuum specific impulse are 
shown in Figs. 14c, 14d, and 14e, respectively. The 
calculated mass flow rates shown in Fig. 14c for the 
3D analyses are lower that the values calculated in 
the 2D axisymmetric analyses. This is because they 
were conducted at a later burn time that corresponds 
to a lower chamber pressure of 4.25 MPa (616 psia). 


The 2D-axisymmetric analyses were conducted at a 
higher chamber pressure of 6.25 MPa (906 psia). 
Thus, the calculated values for the thrust shown 
(Fig. 14d) have similar behavior. The ratio of the 
above is shown as a specific impulse and shown in 
Fig. 14e. The 3D results show lower specific 
impulse because of 3D vectoring of the nozzle. The 
FSM has a higher specific impulse because of the 
higher mass flow rate. 

Summary and Conclusions 

Based on the 2D-axisymmetric analyses 
(Simulation I), the following conclusions have been 
reached : 

• Significant effort has been expended in 
gridding and calculation of the wall y+ in 
conjunction of consistent turbulence models. 

• Convective heat transfer coefficients have 
been calculated using two methods and 4 turbulence 
models. They have been compared with 
measurements, turbulent boundary layer (TBL) 
theory code. Based on the present results of heat 
transfer coefficients and economy, turbulence 
Model 2 on Grid 5 (k-£ RNG with standard wall 
functions as given in Tables 1 and 2) gives the best 
results in relative comparison with the measured 
data [48]. Models 3 and 4 (Table 2) also give 
similar results but time and cost prohibitive. 

Based on the 3D four-degree vectored nozzles 
(Simulation 2), the following conclusions have been 
reached: 

• Chamber: Calculated chamber velocity 

magnitude in the FSM is 25% higher in comparison 
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with the RSRM. This is mainly due to the slight 
opening of the nozzle throat. A side from issues 
concerning erosive burning and maximum chamber 
head-end pressure, the overall behavior of the two 
motors is similar in terms of Mach number and 
static pressure drop. 

Submerged Region: Static pressure varies axially, 
radially, and circumferentially in the submerged 
region of the FSM and the RSRM. A maximum 
circumferential pressure differential of 0.06 atm 
(0.87 psid) and 0.11 atm (1.68 psid) have occurred 
in the RSRM and FSM. respectively. The 
submerged region in the FSM would see a similar 
convective thermal environment in comparison with 
the RSRM. A1 2 Oj concentration and rate of 
accretion are similar in the FSM and the RSRM. 
Nozzle: Slight decrease in the maximum Mach 
number at the nozzle exit. The nozzle in the FSM 
would see a 12% increase in the convective thermal 
environment in comparison with the RSRM. Al 2 Oi 
concentration is much higher in the FSM than in the 
RSRM. AhO< rate of accretion flux in the FSM is 
about 20% higher than in the RSRM. 
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Table 1: Summary of geometry, operating chamber conditions, simulations conducted, grid types and sizes, 
turbulence models used and some pertinent results 


Item 

Simulation 1: RSRM, 2-D Axisymmetric, 
Initial Burn 

Simulation 2a: RSRM, 3D, 
four-degree Vectored Nozzle 

Simulation 2b: FSM, 3D, 
four-deeree Vectored Nozzle 


1.558,61.35 

3.635. 143.11 

3.594. 141.50 


1.37,53.86 

1.364. 53.684 

1.502. 59.146 


3.801, 149.64 

3.792. 149.29 

3.865. 152.173 


1.298. 7.72 

7.10. 7.73 

5.73.6.619 

P» (MPa. psia) 

6.25, 906 

4.25.616 

4.25.616 

r„(K, R) 

3416.04,6148.87 

3368.31,6062.96 

3368.31,6062.96 

Grid 1 

Model 1: k-£, standard, standard wall 
functions, Nx = 1200, Ny = 50, plus sub- 
domains yielding 65,449cells (S, Q), 

(1) AD 

(2) UST = 2644 K 

(3) SST 

1 < y + < 1400, used in [29] 

Model 1: k-e, standard, 
standard wall functions. 

Block 1 (chamber - inner): 
41x57x33, Block 2 (chamber 
-outer): 85x57x33, Block 3 
(nozzle): 101x33x33 yielding 
324,608 cells (S, H) 

(1) SST 
1 <v + < 1800, 

Model 1 : k-e, standard, 
standard wall functions 

Block 1 (chamber - inner): 
41x57x33, Block 2 (chamber - 
outre): 85x57x33, Block 3 
(nozzle): 101x33x33 yielding 
324.608 cells (S, H), SST) 

1 < / < 1 800. 

Grid 2 

Model 1: k-e, standard, standard wall 
functions, (a) Nx = 1200, Ny = 100, plus 
sub-domains yielding 132,323 ceils (S, Q) 

(1) AD 

(2) UST = 2644 K 

(3) UST = 3033 K 

(4) S, Q, SST 

(b) 132000 — > 327599 (S, H, A) 

1 < v + < 700 

Model 1: k-e, standard, 
standard wall functions 

324.608 — » 953,138 (S, H, A) 
( 1 ) SST 

/ < 303 

Model 1 : k-e* standard, 
standard wall functions 

324,608 -> 1,101,636 (S, H, A) 
(1) SST 

v + < 341 

Grid 3 

Model 1: k-e, standard, standard wall 
functions, Nx = 500, Ny = 200, plus sub- 
domains yielding 31370 cells - adapted to 
105161 (S.T, A) 

SST 

1 <v + < 1800 

y* < 2 would require millions 
of cells 

Model 1 : k-e, standard, 
standard wall functions, to 

check exit Mach number 
1,129104 (S, H) 

(1) SST 
v + < 2000 

Grid 4 

71376 (S.T), SST 

9 < y + < 99 >9 



Grid 5 

Model 1: k-e, standard, standard wall 
functions, Nx = 460, Ny = 200, plus sub- 
domains yielding 95550 cells (S, Q) 

SST 

22 < v + < 1 00 



Model 2: k-e, RNG, standard wall 
functions, Nx = 460, Ny = 200, plus sub- 
domains yielding 95550 cells (S, Q) 

(1) SST 
22 < v + < 100 



Grid 6 

Model 3: k-e, RNG, two-layers zonal 
model, 601450 (S, Q) 

(1) SST 
v + < 5 



Model 4: k-e, Standard, two-layers zonal 
model, 601450 (S, Q) 

(1) SST 
v + < 7.5 




A: with grid adaption (shortcoming: calculated \ r and h/Cp profiles have no distinct trend or profile as shown in Figs. 3a and 4b, 
respectively), AD: Adiabatic wall, H: Hexahedral (3D), Q: Quadrilateral (2D), RNG: renormalization group theory, S : structured 
(not triangles), SST : Specified Surface Temperature as shown in Fig. 4a, T: Triangular, UST: Uniform Surface temperature. 













Table 2: Turbulence models used in this study 


Model, Grid 

Viscous 

Model 

K-e Model 

Near Wall Treatment 

Solver 

Calculated wall 

+ 

y 

Model 1, Grid 2 a 

K-e 

Standard 

Standard Wall Functions 

Segregated 

2 < y + < 700 

Model 1, Grid 5 

K-e 

Standard 

Standard Wall Functions 

Segregated 

22 < y + < 100 

Model 2, Grid 5 

K-e 

RNG 

Standard Wall Functions 

Segregated 

22 < y* < 100 

Model 3, Grid 6 

K-e 

RNG 

Two-Layer Zonal Model 

Segregated 

A 

00 

Model 4, Grid 6 

K-e 

Standard 

Two-Layer Zonal Model 

Segregated 

OC 

V 

Model 1, Grid 5 

K-e 

Standard 

Standard Wall Functions 

Coupled 

20 < y + < 100 


d : undersigned for y + ., RNG: renormalization group theory 





Fig. la Local Mach number distribution in the RSRM cavity of the RSRM at ignition time calculated using 

at ignition time, Grid 5 (95,550 quadrilateral cells). turbulence Model 2. 



Fig. lb(l) Local Mach number distribution in the 
cavity of the RSRM at ignition time calculated using 
turbulence Model 1, undersigned fory + . 
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Fig. lb(2) Local Mach number distribution in the 
cavity of the RSRM at ignition time calculated using 
turbulence Model 1. 
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Fig. lb(4) Local Mach number distribution in the 
cavity of the RSRM at ignition time calculated using 
turbulence Model 3. 
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Fig. lb(5) Local Mach number distribution in the 
cavity of the RSRM at ignition time calculated using 
turbulence Model 4. 
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Fig. lc Local static temperature distribution in the 
RSRM at ignition time Grid 5(95,550 quadrilateral 
ceils). 
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Fig. 2a Local static pressure distribution in the RSRM 
at ignition time, Grid 5 (95,550 quadrilateral cells). 
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Fig. 2b Local axial static pressure along the centerline 
of the RSRM at ignition time. 



Fig. 3a Local wall y* along the converging-diverging 
section of in the RSRM nozzle at ignition time. 



Fig. 3b Local walljr* along the converging-diverging 
section of in the RSRM nozzle at ignition time from 
turbulence models 3 and 4 (magnified from previous 
figure). 
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Fig. 3c Comparison of the velocity law of the wall 
(Spalding) with the RSRM nozzle calculations at exit 
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Fig. 4a Local specified surface temperature and 
calculated temperatures in the RSRM at ignition time. 
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Fig. 3d Comparison of the temperature law of the wall Fig. 4b(l) Local convective heat transfer on the nozzle 

with the RSRM nozzle calculations at exit. °f the RSRM at ignition time (Method 1). 
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Fig. 4b(2) Local convective heat transfer on the nozzle 
of the RSRM at ignition time (Method 1, magnified 
from previous figure). 
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Fig. 4b(3) Local convective heat transfer on the nozzle 
of the RSRM at ignition time (Method 2). 
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Fig. 4b(4) Local convective heat transfer on the nozzle 
of the RSRM at ignition time (Method 2, magnified 
from previous figure). 
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Fig. 4c Shear stress profile along the nozzle wall of the 
RSRM at ignition time. 
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Fig, 4d Skin friction coefficient along the nozzle 
centerline of the RSRM at ignition time. 



Fig. 4e Local Reynolds number along the nozzle 
centerline of the RSRM at ignition time. 



Fig. 4f Non-dimensional pressure gradient (P/ + ) along 
the nozzle of the RSRM at ignition time based on 
internal boundary layer. 



Fig. 4g Effective thickness of the viscous sublayer 
(A/) along the nozzle centerline of the RSRM at 
ignition time based on internal boundary layer. 
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Fig. 4h Non-dimensional pressure gradient (/V) along 
the nozzle of the RSRM at ignition time based on 
external boundary layer. 
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Fig. 4j An acceleration parameter (A) along the nozzle 
centerline of the RSRM at ignition time. 
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Fig. 4i Effective thickness of the viscous sublayer (A 2 + ) 
along the nozzle centerline of the RSRM at ignition 
time based on external boundary layer. 



Fig. 5a Rosin-Rammler curve for particle size data. 
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Fig. 5b Local discrete phase concentration on the 
submerged wall of the RSRM at ignition. 
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Fig. 5d Local discrete phase rate of accretion flux in 
the submerged wall of the RSRM at ignition. 
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Fig. 5c Local discrete phase concentration on the Fig* 5e Local discrete phase rate of accretion flux on 

converging-diverging part of the nozzle wall of the the converging-diverging part of the nozzle wall of the 

RSRM at ignition. RSRM at ignition. 
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Fig. 6a Three-dimensional grid (Grid 1 : 325000 Fig. 6d Details of the symmetry plane of the three- 

hexahedral cells) for the RSRM. dimensional grid (Grid 1: 325,000 hexahedrai cells) for 


the FSM. 



Fig. 6b Three-dimensional grid (Grid 1: 325000 Fig. 7a Local flow field in the RSRM on Grid 1. 

hexahedrai cells) for the FSM. 



Fig, 6c Details of the symmetry plane of the three- Fig. 7b Local flow field in the FSM on Grid 1. 

dimensional grid (Grid 1: 325,000 hexahedrai ceils) for 
the RSRM. 
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Fig. 7c Local velocity magnitude distribution at nozzle Fig. 8b Local Mach number distribution in the FSM 

upstream in the RSRM on Grid 1. on Grid 1. 



Fig. 7d Local velocity magnitude distribution at nozzle Fig. 9a Local static pressure distribution in the RSRM 

upstream in the FSM on Grid 1. on Grid 1. 



Fig. 8a Local Mach number distribution in the RSRM Fig. 9b Local static pressure distribution in the FSM 

on Grid 1. on Grid 1. 
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Fig. 9c Local static pressure gradient distribution (five 
planes in the submerged region and parallel to inlet) in 
the RSRM on Grid 1. 
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Fig. 10a Local values of wall y + in the RSRM on Grid 
1. 
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Fig. 9d Local static pressure gradient distribution 
(four planes in the submerged region and parallel to 
inlet) in the FSM on Grid 1. 

O RSRM, 3D, 4 deg □ FSM, 3D, 4 deg Os from [7] 
x 9 s from [7] X 19 s from [7] • 114 s from [7] 



Fig. 9e Local static pressure gradient in the RSRM 
and FSM on Grid 1 along with cold flow experiments 
of Whitesides et al. [7] 


Fig. 10b Local values of wall y in the FSM on Grid 1 















Fig. 11a Local convective heat transfer in the RSRM 
on Grid 1. 


Fig. 12a Local discrete phase concentration in the 
RSRM on Grid 1. 
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Fig. lib Local convective heat transfer in the FSM on 
Grid 1. 


Fig. 12b Local discrete phase concentration in the 
FSM on Grid 1. 
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Fig, 13a Local discrete phase rate of accretion flux in 
the RSRM on Grid 1. 
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Fig. 13b Local discrete phase rate of accretion flux in 
the FSM on Grid 1. 
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Fig. 14a Chamber pressures. 



Fig. 14b Nozzles expansion area ratio. 



Fig. 14c Calculated mass flow rates. 
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Fig. 14d Calculated vacuum thrust. 
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Fig* 14e Calculated vacuum specific impulse. 
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